Population structure and demographic history of the gastropod Thaisella chocolata (Duclos, 1832) from the Southeast Pacific inferred from mitochondrial DNA analyses

Abstract The present‐day population structure of a species reflects the combination of oceanographic currents, life‐history traits, and historical events. However, little is known about the mechanisms that have shaped the gene lineage distribution of marine species inhabiting the Southeast Pacific. Here, we provide a comprehensive phylogeographical study of a species distributed along the Southeast Pacific coastal region by analyzing the endemic gastropod Thaisella chocolata (Duclos, 1832). Sequencing of mitochondrial cytochrome c oxidase subunit 1 (CO1) and 16S rRNA revealed strikingly high haplotypic nucleotide and genetic diversity but a lack of significant population differentiation within the survey area. In addition, a star‐shaped phylogeny and significantly negative Tajima's D and Fu's Fs tests of neutrality suggested historical occurrence of rapid demographic expansion. Mismatch distributions and Bayesian inference analyses also confirmed T. chocolata to have undergone two ancestral demographic expansions. Calculations suggested that these expansions began in the lower and middle Pleistocene epoch, likely due to continental shelf development and climatic conditions. These findings could help establish a genetic baseline for T. chocolata as the first step toward sustainable spatial management of this species, as well as understand this species’ response to future climate change.

graphical study of a species distributed along the Southeast Pacific coastal region by analyzing the endemic gastropod Thaisella chocolata (Duclos, 1832). Sequencing of mitochondrial cytochrome c oxidase subunit 1 (CO1) and 16S rRNA revealed strikingly high haplotypic nucleotide and genetic diversity but a lack of significant population differentiation within the survey area. In addition, a star-shaped phylogeny and significantly negative Tajima's D and Fu's Fs tests of neutrality suggested historical occurrence of rapid demographic expansion. Mismatch distributions and Bayesian inference analyses also confirmed T. chocolata to have undergone two ancestral demographic expansions. Calculations suggested that these expansions began in the lower and middle Pleistocene epoch, likely due to continental shelf development and climatic conditions. These findings could help establish a genetic baseline for T. chocolata as the first step toward sustainable spatial management of this species, as well as understand this species' response to future climate change.

K E Y W O R D S
16S rRNA, cytochrome oxidase subunit 1, genetic diversity, phylogeography, population structure, Thaisella chocolata

T A X O N O M Y C L A S S I F I C A T I O N
Life history ecology

| INTRODUC TI ON
An essential component of long-term fisheries management is the determination of a species' genetic population structure, which commonly reflects the combination of oceanographic currents, life-history traits, and historical events . In the Southeast Pacific, climatic and geological events especially have played significant roles in modifying species assemblages and species ranges Cárdenas, Silva, et al., 2009). However, distinct from the Atlantic and the Northeast Pacific, how oceanographic factors and historical events have helped shape the region's rich biodiversity is poorly documented. One of the most important geological events in the Southeast Pacific is the appearance of the modern Humboldt Current System (HCS), one of the four major global eastern boundary currents (Hill et al., 1998;Zinsmeister, 1978). The HCS is a strong coastal upwelling ecosystem that derives from the West Wind Drift (WWD) at around 40° S and flows northwards along the Pacific eastern seaboard as part of the Coriolis-force-induced South Pacific gyre (Gutiérrez et al., 2016).
The HCS comprises four major currents: the combination of the Peru and Chile coastal currents that exist over the shelf area and reach only to shallow depths; the Peru-Chile current found at the boundary between the upwelling and downwelling regions; the Peru-Chile Countercurrent, a poleward flow that brings cold and nutrient-rich water to the surface from depths below 150 m; and the Peru-Chile undercurrent typically found at depths between 50 and 4400 m (Karstensen & Ulloa, 2009). By creating regional geographic barriers, these currents can influence local adaptation in marine species populations; for example,  proposed that the two spatial genetic patterns found in the marine gastropod Concholepas concholepas along the length of the HCS could be an effect of the WWD current, which splits at approximately 43° S into the equatorward Peru Current and the pole-ward Cape Horn Current. In addition, being geographically situated near the eastern terminus of the Pacific equatorial wave guide, the HCS bears the immediate brunt of the El Niño-Southern Oscillation (ENSO), an interannual climate variability that causes strong fluctuations in species abundance and dynamics (De Oliveira et al., 2009;Martínez et al., 2003). During the warm phase of an ENSO event, the wind field changes at the equator and the warm waters of the western equatorial Pacific region propagate to the eastern Pacific by means of Kelvin waves that radiate poleward along the eastern boundary in coastal trapped waves (Philander, 1983;Thiel et al., 2007). The HCS is thereby subjected to intermittent drastic environmental rearrangements more extreme than those known to occur in any other sizeable marine ecosystem on Earth (Bakun & Weeks, 2008).
Because the HCS has been a permanent oceanographic feature in the Southeast Pacific throughout most of the Cenozoic, the geologic record along the coasts of Peru and Chile should contain a detailed history of past climatic events (Zinsmeister, 1978). In fact, the marine species inhabiting this region evidence genetic signatures that reflect demographic changes driven by geological and oceanographic changes during the Tertiary and early Quaternary Cárdenas, Silva, et al., 2009). For example, prior to the development of the West Antarctic ice sheet, the prevailing westerly winds would have been far weaker, with the net result of a southward shift of the Polar Front and a corresponding southward displacement of the WWD (Zinsmeister, 1978). This more southerly placement may have led to prevailing wind patterns along the Southeast Pacific coasts to promote the development of a strong, warm, south-flowing, ENSO-like countercurrent along Peru and Chile (Zinsmeister, 1978). The presence of subtropical molluscan faunas as far south as the Golfo de Penas, Chile, in the Middle and Early Miocene, and the presence of many extralimital species in the Antofagasta-Mejillones area, Chile, during the Marine Isotope Stage 11 interglacial (MIS-11) suggests that ENSO conditions indeed existed (Berger & Wefer, 2003;Zinsmeister, 1978). It is envisioned that during the MIS-11 interglacial (426-396 ka before present [BP]), the climatic regime of northern Chile was characterized by significantly warmer air temperature, increased solar radiation, reduced cloud cover, and enhanced strong, quasi-periodical, ENSO-like conditions (Ortlieb et al., 2003;Tzedakis et al., 2022). If this hypothesis is correct, warmer conditions would have intermittently existed much farther south along South America (Zinsmeister, 1978). As mentioned above, ENSO events may become biological disasters affecting populations' genetic signatures (De Oliveira et al., 2009;Martínez et al., 2003). If a single event can have such an impact, the question arises as to what the cumulative effects may be over a longer time scale (Ibáñez et al., 2012).
Glacial and interglacial periods and consequent eustatic sea-level oscillations during the Late Pleistocene also influenced the distribution and population sizes of temperate marine species inhabiting the Southeast Pacific . The interval spanning from the Last Glacial Maximum (LGM, 21-18 ka BP) to the preindustrial era represents the most recent large-scale reorganization of the climate system (Osman et al., 2021). During this time, an abrupt decline in sea level, to about 120-140 m below current sea level, induced a massive exposure of the continental shelf that left genetic signatures at population and species levels congruent with population decline (Lambeck et al., 2002;Salvatteci, 2013). In the Southeast Pacific, for example, analysis of laminated sediments of the Peruvian continental shelf revealed minimum abundances of oceanic species during the LGM (Salvatteci, 2013). Then, at the beginning of the Younger Dryas period (12.9-11.8 ka), during which the sea level rose, the primary productivity off Peru increased (Rein et al., 2005). These fluctuating paleoclimatic conditions can be expected to have consequences on both demography and population structure Echevin et al., 2011). For the Southeast Pacific, however, most biogeographical studies to date have focused on the description of species population structures without a deep phylogeographical assessment of the mechanisms that could have generated the observed patterns. This limits our understanding of the historical and contemporary processes that led to the present-day distribution of the species in question. To address this issue, in this study, we selected a species with a population history that can be traced back to the latest part of the Pliocene epoch and that displays a geographical range covering the entire Peruvian coastal region, including the transition zone between the Peruvian and the Panamanian biogeographical provinces: the marine gastropod Thaisella chocolata (Duclos, 1832) (DeVries, 2007Ganoza-Chozo et al., 2021;Osorio-Ruíz, 2002;Tarazona et al., 2003). On account of being endemic to the Southeast Pacific coast and having a life history that includes a planktonic larval stage of 4 months (Ganoza-Chozo et al., 2021;Osorio-Ruíz, 2002;Romero et al., 2004), T. chocolata, an inhabitant of the rocky intertidal zone of South America, is ideal for studying the historical and contemporary process that led to its present-day distribution and the effects of sharp contrasts among biogeographical areas. Moreover, T. chocolata is a heavily exploited gastropod, with extraction records going back to the 1970s (Flores et al., 1998;Ganoza-Chozo et al., 2021;Promperú, 2019). Therefore, identification of potential subpopulations of this species, as in this study, may provide an important regional context in which to frame marine management or conservation efforts. Conservation efforts are especially necessary since the Peruvian population of T. chocolata has rapidly decreased due to fishing pressure and anthropogenic disturbances (Alfaro-Mudarra, 2020).
Even though marine populations may consist of ample numbers of individuals to start showing declines in genetic diversity measurable with feasible sample sizes (Ryman et al., 1995), the effective population size, which determines the genetic properties of a population, may be orders of magnitude smaller than the census population size, for example, millions of individuals may be equivalent to an effective population size of only hundreds or thousands (Hare et al., 2011).
The notion that marine species' populations lose genetic diversity is thus possible despite their large spawning stock biomasses, which could potentially result in reduced adaptability and productivity (Hauser et al., 2002). Considering that the Peruvian T. chocolata fishery is based on stocks that have already shown a reduction in population size and that the average catch for this species in the last decade is around 28% less compared with the average record of marine fishery landings reported since 1975 (Alfaro-Mudarra, 2020;Flores et al., 1998;PRODUCE, 2021), fishery may have already caused considerable changes in genetic diversity. Determining the population's genetic diversity is, therefore, paramount to conservation and fishery decision-making.
Genetic tools are a powerful means of examining demographic history and population dynamics; in particular, mitochondrial DNA (mtDNA) has an effective population size of approximately onequarter that of nuclear markers, which allows the patterns of recent historical events to be recovered without extensive sequencing efforts (Hurst & Jiggins, 2005). Moreover, although doubly uniparental inheritance (DUI) of mtDNA occurs in some invertebrate species (Ladoukakis & Zouros, 2001;Soroka, 2020), there is currently no evidence of DUI in gastropods (Gusman et al., 2017). Thus, maternal inheritance and lack of recombination also make mtDNA an appropriate marker for tracing recent evolutionary history, including colonization events and population bottlenecks, without the confounding effects of biparental inheritance and recombination inherent in nuclear DNA (Avise, 2000;Harrison, 1989;Pakendorf & Stoneking, 2005). The use of mtDNA in the study of demographic history has been especially revealing. For example, sequencing of cytochrome c oxidase subunit 1 (CO1) and the mitochondrial control region has allowed the effects of glacial and interglacial episodes to be determined among species of global commercial importance in the Southeast Pacific region, such as Dosidicus gigas (Ibáñez & Poulin, 2014), Scomber japonicus , and Trachurus murphyi Cárdenas, Silva, et al., 2009). Given a clear understanding of pelagic populations' life histories and dynamics, such analyses can yield helpful information for understanding the species' response to future climate change. In addition, fishery management initiatives provided with such understanding are much more likely to be able to protect aspects of biocomplexity such as spawning areas and local adaptations.
In the present study, we present the phylogeographical pattern of the marine gastropod T. chocolata within its native range as determined from the analysis of two mitochondrial markers, cytochrome oxidase subunit 1 (CO1) and 16S rRNA. Sample collection covered most of the present-day distributional range of the species, with populations representing each of the three biogeographical regions described for the Peruvian coast: northern, central, and southern.
Using population genetic analyses, we first utilized haplotype network computations and molecular-based methods to determine the gene lineage distribution of T. chocolata along its range. Second, we aimed to decipher the biogeographic pathways that T. chocolata experienced, considering the geological and climatic history of the region. To achieve these aims, our study combined coalescent theory, multivariate analyses, Bayesian clustering, and Bayesian skyline reconstruction to facilitate a more realistic understanding of this gastropod's demographic history. Per-location, sample sizes ranged from 18 to 22. For each captured snail, a tissue sample was collected from the foot muscle. Total genomic DNA was extracted from samples by means of proteinase K digestion followed by cetyltrimethylammonium bromide (CTAB) extraction according to a protocol adapted from Zuccarello and Lokhorst (2005). DNA quality and quantity were checked using the Epoch microplate spectrophotometer (BioTek Instruments); afterward, the concentration was adjusted to 50 ng μl −1 . Partial regions of the mitochondrial CO1 gene and 16S rRNA were amplified using universal primers, HCO2198 and LCO1490 for CO1 (Folmer et al., 1994) and 16Sar and 16Sbr for 16S rRNA (Palumbi, 2005). All novel CO1 and 16S rRNA sequences generated in this study have been submitted to the GenBank database (www.ncbi.nlm.nih.gov/ genbank) under accession numbers OK087138-OK087293 and OK094932-OK095087 (Tables S1 and S2).

| Genetic diversity and population structure
The nucleotide sequences obtained for T. chocolata mitochondrial CO1 and 16S rRNA were aligned in MAFFT 7 (Katoh & Standley, 2013) using default settings and trimmed to the shortest sequence. The number of haplotypes (N H ), haplotype diversity (h), number of segregating sites (S), and nucleotide diversity (π) among samples were calculated using DnaSP 6.12.03 (Rozas et al., 2017). The number of private haplotypes unique to each population (N P ) was determined according to the haplotype list for each population generated in Arlequin 3.5.2.2 (Excoffier & Lischer, 2010). Haplotype networks were constructed with the TCS statistical parsimony algorithm implemented in PopART 1.7 (Leigh & Bryant, 2015). To assess population structure, we estimated Wright's F-statistic (F ST ) (Weir & Cockerham, 1984) using Arlequin with statistical significance determined on the basis of 10,000 permutations. PhiST (Φ ST ), an estimator of genetic differentiation independent of mutation rate (Kronholm et al., 2010), was also computed using the R package diveRsity (Keenan  et al., 2013), and statistical significance examined using 10,000 bootstrap iterations. Hierarchical genetic population structure was tested for with the analysis of molecular variation (AMOVA) implemented in Arlequin. To further evaluate and visualize that population genetic structure with respect to geography, principal coordinates analysis (PCoA) was conducted using DARwin 6.0.21 (Pirrier & Jacquemoud, 2006).

| Haplotype accumulation curves
A nonparametric sample size estimate was determined using the R package spider (Brown et al., 2012). We employed random permutation subsampling via the haploAccum() function and by using the Chao 1 estimator (Chao, 1984) with the chaoHaplo() function. The latter estimates the minimum sample size capable of accounting for all haplotype diversity in a genomic dataset based on both the number of haplotypes present and the numbers of singleton and doubleton sequences (i.e., those occurring once and those appearing twice) (Brown et al., 2012;Robalo et al., 2020).

| Population demographic history
The population demographic history and evolutionary neutrality of T. chocolata were assessed using two tests, Fu's Fs (Fu, 1997) and Tajima's D (Tajima, 1989), both implemented in Arlequin. p-values were calculated by means of 10,000 coalescent-based simulations.
A mismatch distribution analysis using a sudden expansion model in DnaSP was conducted to compare the distribution of pairwise differences between haplotypes (Harpending, 1994;Rogers & Harpending, 1992). A Bayesian skyline plot (BSP) was also employed to examine historical demographic fluctuations since the time of the most recent common ancestor; this plot was created using Beast 3.6.4 , hosted on the CIPRES Science Gateway (Miller et al., 2010). For both regions of interest, the BSP was built using a coalescent Bayesian skyline model under an HKY substitution model with empirical base frequencies and a strict molecular clock; however, certain parameters differed. For CO1, the substitution rate was fixed to a mutation rate of 1.57% My −1 previously applied in evolutionary studies of invertebrates (Wilke et al., 2009), and in the coupled MCMC analysis, the delta temperature was set to 0.1. The analysis was run in triplicate for 700 million generations with a sampling frequency of 70,000 and burn-in of 10%. For 16S rRNA, the substitution rate was fixed to 1.4% My −1 , a value obtained from the combined CO1 and 16S rRNA data (Morrison et al., 2004). The analysis was run in triplicate for 500 million generations with a sampling frequency of 50,000 and burn-in of 10%. In both cases, the triplicate analyses were combined using LogCombiner . Convergence of MCMC chains was assessed by visual examination in Tracer (Rambaut et al., 2014) of the log trace of each posterior distribution showing an ESS greater than 200. Files from all BSP analyses were collapsed into graphs using Tracer and R (R Core Team, 2020).

| Population genetic analyses
The mitochondrial gene CO1 (508 bp long) was amplified and sequenced from 156 individuals of T. chocolata collected along the northern, central, and southern regions of the Peruvian coast.
Across the total sample, the number of haplotypes was 129, while the number of polymorphic sites was 107, of which 37 were singleton variable sites and 70 parsimoniously informative sites.
This remarkable mtDNA diversity was further reflected by a haplotype diversity value close to its maximum of 1, which indicates that in the overall population, the probability of two individuals from the same locality sharing the same haplotype is <0.001%.
Population-wise, haplotype and nucleotide diversities were also very high (h = 0.963 ± 0.033 to 1.000 ± 0.019; π = 0.013 ± 0.001 to 0.016 ± 0.001) ( Table 1). Even though almost every haplotype was a singleton, pairwise comparisons revealed that one of the haplotypes (Hap8) was present throughout the sampling time. Mitochondrial 16S rRNA (410 bp long) was also successfully amplified and sequenced from 155 individuals of T. chocolata. The number of haplotypes detected across the sample was 33, while the number of polymorphic sites was 38, of which 28 were singleton variable sites and ten parsimoniously informative sites. All populations shared at least one haplotype, and the most common haplotype (By17) was present throughout the sampling time. Population-wise, haplotype and nucleotide diversities for this locus ranged from 0.284 ± 0.128 to 0.763 ± 0.103 and from 0.0007 ± 0.0004 to 0.0039 ± 0.0009, respectively (Table 1).

| Population structure and phylogenetic analyses
We next investigated whether population genetic structure through time and space could contribute to the mtDNA hyperdiversity observed in T. chocolata. At first glance, the high number of private haplotypes (Table 1)  rRNA genomic datasets did not reveal any visually discernible patterns (Figure 4b,d). Instead, all populations appeared to almost completely overlap with each other. Finally, evaluation of the population structure of the native range collections via AMOVA confirmed an absence of genetic structure between populations; differences between areas inferred from the CO1 and 16S rRNA genomic datasets accounted for 0.86% y 0% of the total variance, while 99.14% and 100% was attributable to differences among individuals within each location (Table 2). These results indicate panmixia of T. chocolata across the geographical study area.

| Haplotype accumulation curves
The haplotype accumulation curves for both CO1 and 16S rRNA genomic datasets failed to reach the asymptote (Figure 5a,c), indicating that only part of the actual genetic diversity was captured.
Estimates of the sample sizes needed to capture the total estimated TA B L E 1 Summary statistics for mitochondrial markers cytochrome c oxidase subunit 1 (CO1) and 16S rRNA in Thaisella chocolata

| Demographic analysis
Fu's Fs and Tajima's D neutrality tests were performed for the CO1 and 16S rRNA haplotypes (Table 3)

| DISCUSS ION
Fishery management to date has largely been concerned with the immediate resource of interest-the abundance and size of the fish available for harvesting (Allendorf et al., 1987). Little attention has been directed toward understanding the genetics of their populations, despite such understanding being crucial for determining the patterns of genetic connectivity that ultimately define the geographic structure of a species (Allendorf et al., 1987;Thiel et al., 2007). In the present study, the marine snail T. chocolata was found to exhibit remarkably high haplotype and nucleotide genetic diversity throughout its Peruvian geographical range, Animal mitochondria have been shown to be inefficient, if not outright lacking, in repairing various types of DNA damage, a factor that by itself could contribute significantly to the high mutation rate observed for mtDNA (Brown et al., 1979;Clayton et al., 1974;Palit & Ngili, 2018). Thomas et al. (2010) also showed that in invertebrate eumetazoan species with shorter generation times, protein-coding F I G U R E 2 Haplotype network and phylogenetic tree based on the cytochrome c oxidase subunit 1 (CO1) genomic dataset. (a) The unrooted TCS haplotype network was constructed using 129 CO1 haplotypes from Thaisella chocolata. Each line between points represents a single mutational step. Circle size indicates haplotype frequency. TCS network was generated using PopART. (b) Phylogenetic tree reconstructed using 156 CO1 sequences. Bayesian inference analysis was performed for 100 million generations using MrBayes. Stramonita biserialis was used as an outgroup.

F I G U R E 3
Haplotype network and phylogenetic tree based on the 16S rRNA genomic dataset. (a) The unrooted TCS haplotype network was constructed using 33 16S rRNA haplotypes from Thaisella chocolata. Each line between points represents a single mutational step. Circle size indicates haplotype frequency. TCS network was generated using PopART. (b) Phylogenetic tree reconstructed using 155 16S rRNA sequences. Bayesian inference analysis was performed for 100 million generations using MrBayes. Stramonita biserialis was used as an outgroup. vertebrates (Bazin et al., 2006), the pattern of low genetic diversity in mtDNA does not apply to all species. Indeed, mtDNA hyperdiversity has been reported in many temperate gastropod species and TA B L E 2 Hierarchical analysis of molecular variance (AMOVA) for the cytochrome c oxidase subunit 1 (CO1) and 16S rRNA genomic datasets  highly diverse planktonic-dispersing marine invertebrates (Fourdrilis et al., 2016). In the marine gastropod Melarhaphe neritoides, especially, Fourdrilis et al. (2016) reported the gene CO1 to have a mutation rate of as high as 1.99 × 10 −4 per site per generation, a value 1000-to 10,000-fold higher than commonly estimated for mtDNA in metazoans from other phyla. This suggests that mtDNA hyperdiversity may be more common across eukaryotes than currently known. Thus, the high mutation rate of mtDNA could be the main factor driving the observed CO1 hyperdiversity in T. chocolata. It is likely, however, that the rate of variation in mtDNA evolution is attributable to a more complex combination of factors than a high mutation rate alone. A second hypothesis suggests that the observed mitochondrial hyperdiversity could be explained by adaptive evolution. By evaluating the ratio of nonsynonymous to synonymous changes within and between species, Bazin et al. (2006) determined invertebrate mtDNA loci to exhibit a significant shift toward values less than one, consistent with adaptive evolution. In M. neritoides, Fourdrilis et al. (2018) found that CO1 is under purifying selection and that among mitochondrial genes, it is most strongly subject to selection. Given this non-neutrality, mtDNA diversity would, in many instances, reflect the time since the last occurrence of a selective sweep rather than population history and demography (Bazin et al., 2006). We should mention, however, that Quesada et al. (1998) pointed out that an excess of replacement substitutions could be ex-  Cárdenas, Silva, et al., 2009;Fernández-Iriarte et al., 2020;Ibáñez et al., 2011;Ibáñez & Poulin, 2014). The absence of geographic genetic structure observed in these species seems to obey the apparent lack of major oceanographic barriers within this region. It is also feasible that the combination of northward (Peru current) and southward (Peru-Chile countercurrent) oceanographic flows allow effective population mixing. In addition, oceanographic features of the HCS, such as the ENSO, could also contribute to genetic homogeneity . During the El Niño warm phase of the ENSO pattern, warm water masses moving southward cause the transport of larvae to the south, whereas, during the cool La Niña phase, retreat of surface equatorial waters toward the north causes larval transport northward (Días & Ortlieb, 1993;Paredes et al., 1998). Together, the constant northward transport of the Humboldt current and semi-regular southward transport mediated by the ENSO would create a bidirectional gene flow that prevents any degree of genetic differentiation for most species inhabiting the Peruvian coast (Barahona, 2017).

Source of variation
The genetic homogeneity observed in T. chocolata can also be rooted in its high potential for dispersal. The species has an intracapsular development of 49 days and a larval development period of 120 days (Romero et al., 2004); this long spawning time and   (Marko et al., 2010); in species from the northwestern Pacific such as the gastropod Thais clavigera (Guo et al., 2015); and in codistributed rocky-shore species such as Concholepas concholepas Cárdenas, Silva, et al., 2009).
Hence, events predating the LGM seem to have been more critical in shaping the demographic histories of rocky-shore species (Marko et al., 2010).
The first wave of population expansion determined for T. chocolata, at around 800,000-700,000 years BP, also coincides with the most extreme low and high stands of the mid-Pleistocene Revolution (Scourse, 2013). By 800,000 years BP, glacial cycles with sea level amplitudes of more than 100 m were occurring, mostly with periods on the order of 100,000 years (Berger et al., 2016). During these glacial-interglacial cycles, significant eustatic movements took place, including sea-level lowering of at least several tens of meters during cold events and up to 100-140 m during full glacial episodes (Rabassa et al., 2005). These caused the South American coastlines to move westward by 100 km between 7° and9° S and by 20 km at 4°40′ S, resulting in the exposure of a large part of the South Pacific continental shelf (Ortlieb & Machare, 1989). Such substantial marine regressions probably caused habitat expansions for benthic invertebrates such as T. chocolata. In fact, DeVries (2007) showed that the smooth, refractory, chocolate-brown outer shell layer of T. chocolata emerged as coastal water was cooling, submergent shorelines were rising, eustatic sea level was undergoing rapid and marked changes, and Peruvian shorelines were becoming straighter and more exposed. Moreover, Valdovinos et al. (2003) showed that the relative importance of factors controlling species diversity might have shifted in high latitudes along the southeastern Pacific, with available shelf area playing a more prominent role than temperature (Valdovinos et al., 2003). Hence, the first population expan-  (Maggs et al., 2008). Leading-edge populations are those whose demographic and allelic composition stochasticity are predictably stronger (Robalo et al., 2020) and thus are typically affected by founder events. Also, in the wake of a bottleneck or founder event, the geographic space becomes structured into regions of low diversity separated by sharp frequency gradients (Excoffier & Ray, 2008;Ibáñez et al., 2011;Robalo et al., 2020). In our datasets, haplotypes present in the leading-edge patches were not shared with sites at the center of the distribution, and high haplotypic and nucleotide diversity values were observed for all sites.
T. chocolata therefore seems to have retained exceptional levels of genetic diversity through its historical demographic expansions, an event that seems not consistent with a genetic bottleneck or founder effect.  (Cossu et al., 2021;Laikre et al., 2005). The genetic identification of the stock structure, however, is only the first step in implementing a sustainable management plan. Fishery managers should also consider estimating population abundance and population dynamics rates to reach appropriate levels of fishing mortality to maintain population sustainability and productivity (Cochrane, 2002). In fact, the Code of Conduct for Responsible Fisheries specifies that the states should adopt appropriate measures, based on the best scientific evidence available, to maintain or restore stocks at levels capable of producing maximum sustainable yield (FAO, 1995). If experimental data cannot be obtained to reach this objective, species dynamics can be characterized using individual-based simulations such as those implemented by Strand (2002) and already applied to fisheries management (IWC, 2007). These simulations represent natural conditions in an explicitly stochastic fashion by treating individuals as autonomous units that move through the environment, experience demographic transitions, and interact with each other based upon computer-generated random numbers (Strand, 2002). Managers should also consider evaluating the impacts of the changes in any of the ecosystem's biological, chemical, or physical components on the resource population and community (Cochrane, 2002). For example, in some South American countries, the use of butyltin (BT), an organotin compound used in antifouling paints, has not been regulated even though it was banned by the International Maritime Organization because it triggers imposex (IMO, 2000).
Imposex is a type of sexual abnormality wherein male sex organs develop on female gastropods (Zou, 2019). Although BT causes female infertility and inhibits the release of capsules (often containing eggs) (Gibbs & Bryan, 1986;Stroben et al., 1995) Ultimately, management plans should include threshold criteria for detecting biologically significant changes (Flanagan et al., 2018).
A practical approach could consist of set trigger points throughout the range of an indicator variable to ensure that management action is initiated before a crisis point is reached. For example, Flanagan et al. (2018) suggest that a continuous decline in allelic richness at putatively adaptive loci or observation of low survival or fecundity over multiple sampling periods may trigger a management intervention such as genetic rescue to increase allelic richness or fitness. Achieving sustainable management, therefore, will require a proactive approach to seeking ways to optimize the benefits derived from the resource.

| CON CLUS IONS
Measuring dispersal based on genetic patterns is essential to conservation and species management, as dispersal directly determines a species' response to selection and its capacity to respond to disturbances in a changing environment (Bradbury & Bentzen, 2007). Determining dispersal pathways through an understanding of genetic connectivity is, therefore, an essential prerequisite for devising effective fishery management strategies. Based on the analysis of two mitochondrial markers, CO1 and 16S rRNA, the present study examined the genetic connectivity of the gastropod T. chocolata across the Peruvian coastal region.
Our findings support the existence of a genetically homogeneous population within the survey area, which reconciles biological and management criteria. Nevertheless, as population genetic inferences based only on mtDNA can be misleading Cárdenas, Silva, et al., 2009), these results must be treated with caution. For example, besides population expansion, negative values of Tajima's D can reflect other demographic scenarios such as selective sweeps.
In both cases, the lineages in the genealogy are forced to coalesce at the time of the selective sweep or bottleneck-the average number of pairwise differences is decreased compared to the number of segregating sites, leading to negative Tajima's D values (Nielsen, 2001). It is therefore difficult to distinguish both scenarios when a single locus is considered. Using additional markers, a selective sweep in mtDNA can be differentiated from a demographic expansion if a low diversity in the whole genome is revealed (Ibáñez et al., 2011). Hence, further analyses using nuclear DNA markers or whole-genome analysis could provide a complete perspective on population demography for this species. The accumulation of biological and ecological data also remains crucial, as such data are needed to determine the extent of migration, which will help to better define the mode of gene flow between populations (Sekino & Hara, 2001). Demographic history analyses in the present work also revealed that T. chocolata underwent two waves of population expansion, most likely during the lower and middle Pleistocene epoch. According to our data, T. chocolata retained exceptional genetic diversity levels following population expansion.
This finding could provide helpful information for understanding this species' response to future climate change. The analysis of population connectivity in T. chocolata, therefore, can help establish a genetic baseline for this species in order to better maintain the equilibrium of the population while it is subjected to reasonable exploitation.

ACK N OWLED G M ENTS
We thank Gustavo Vasquez, Félix Uchpa, and Roy López for the samples provided for this study.

FU N D I N G I N FO R M ATI O N
This research did not receive any specific grant from the public, commercial, or not-for-profit funding agencies.

CO N FLI C T S O F I NTE R E S T
The authors report no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in the GenBank database (www.ncbi.nlm.nih.gov/genbank) under accession numbers OK087138-OK087293 and OK094932-OK095087. CO1 and 16S rRNA genomic datasets used during this study are also stored and available at https://doi. org/10.5061/dryad.jsxks n0d1.